Universality in three-dimensional Ising spin glasses: A Monte Carlo study 
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We study universality in three-dimensional Ising spin glasses by large-scale Monte Carlo sim- 
ulations of the Edwards-Anderson Ising spin glass for several choices of bond distributions, with 
particular emphasis on Gaussian and bimodal interactions. A finite-size scaling analysis suggests 
that three-dimensional Ising spin glasses obey universality. 
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I. INTRODUCTION 

One of the cornerstones of the theory of critical phe- 
nomena is the concept of universality, according to which 
the values of many quantities, such as critical exponents, 
do not depend on microscopic details but only on a few 
broad features such as the dimensionality of space and 
the symmetry of the order parameter. Universality fol- 
lows from renormalization group (RG) theory, according 
to which many interactions that could be added to the 
Hamiltonian are "irrelevant," i.e., do not change the crit- 
ical behavior. The "e-expansion" implementation of RG 
has been very successful in predicting which perturba- 
tions are actually relevant and which are irrelevant. At 
least for pure systems, and disordered systems without 
frustration, numerical simulations seem to be consistent 
with universality. 

However, for systems with both frustration and disor- 
der, known as spin glasseSfi the situation is less clear. 
On the one hand, e-expansion calculations as well as 
high-temperature series expansions^- for spin glasses im- 
ply universality, and, in the opinion of the authors of 
this paper, there is no a priori reason why universality 
should be less valid for spin glasses than for pure sys- 
tems. However, as we shall see below, numerical results 
so far have not been compelling in favor of universality 
and some groups^i^i^iSii even claim explicitly that univer- 
sality is violated^ Unfortunately, it is difficult to obtain 
accurate critical exponents because there are significant 
corrections to scaling, there are long equilibration times 
in Monte Carlo simulations that limit the available sys- 
tem sizes, and all quantities need to be averaged over 
many realizations of the disorder in order to have small 
enough error bars. 

According to universality, the range of the interac- 
tions is irrelevant, as long as it is finite, and so, for 
example, adding next-nearest-neighbor couplings to a 
nearest-neighbor model will not change the critical be- 
havior. However, random systems are characterized not 
just by the strength of first-neighbor, second-neighbor, 
etc., interactions, but by the distributions of these (ran- 
dom) quantities. Hence, even if one restricts oneself to 
nearest-neighbor interactions, there are many different 
models characterized by different distributions which are 



expected to be in the same universality class. In this pa- 
per we attempt to answer, through careful simulations, 
whether this expectation is true for Ising spin glasses in 
three dimensions. 

Many groups have estimated critical exponents for spin 
glasses for the Edwards- Anderson (EA) Ising spin glasaS 
in three dimensions for different disorder distributions, 
mainly the Gaussian and bimodal (± J) models. In Ta- 
ble m we present a summary of these results. We also 
present results from a recent studylS which uses a three- 
dimensional diluted Ising spin glass (45% bond occupa- 
tion). The advantage of such a model is that, due to the 
dilution, cluster algorithms ^"'^ can be used to study larger 
system sizes and that corrections to scaling seem to be 
small^*' when the bonds are drawn from a bimodal distri- 
bution. The data in Tabled show clearly that there is a 
large spread in the estimates of the different critical ex- 
ponents obtained using several different methods, such as 
series expansions, nonequilibrium relaxation approaches, 
and finite-temperature Monte Carlo methods combined 
with a finite-size scaling (FSS) analysis. The spread in 
the different estimates of the critical exponents and 77 
is easily visualized in Fig. ^ where we plot 77 versus v. 
As mentioned above, some groups claim universality is 
violatedj^tiiSiSii but most of the papers do not make this 
claim, probably because the error bars on the data are 
usually large. In this paper we aim to test universality 
in equilibrium more precisely by reducing the error bars. 

A major problem with reducing error bars in critical 
exponents is the presence of corrections to FSS, which 
means that the scaling expressions used to determine ex- 
ponents do not work well for small system sizes. For 
pure systems, several methods have been proposed in an 
attempt to reduce errors caused by corrections to scaling. 

First, try to eliminate the leading correction. In this 
approach, the model is altered until the operator which 
gives the leading correction does not appear in the Hamil- 
tonian. This means that its effects will not be felt in any 
calculated quantity. For the three-dimensional Ising fer- 
romagnet, using a "soft-spin" model rather than "hard" 
±1 spins, and varying the coefficient of the fourth-order 
term in the Hamiltonian, it was possible to eliminate the 
leading correction to FSS and obtain high-precision val- 
ues for the critical exponentsi^S We have attempted to 
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FIG. 1; (Color online) Graphical representation of different 
estimates of the critical exponents 77 and v (taken from Table 
UJ. The centers of the ellipses represent the different esti- 
mates, and the major axes represent the error bars of the 
estimates. The data show a considerable spread. The thick 
lines represent the estimates from the current study. Note 
that the data from the current study for bimodal (solid lines) 
and Gaussian (dashed lines) distributions agree within error 
bars thus providing evidence for universality. The dotted line 
represents data for a random-anisotropy Heisenberg model in 
the strong-anisotropy limit (Ref. which is expected to be 
in the same universality class as the Ising spin glass. 

do this for the spin glass by (i) choosing different disor- 
der distributions, and (ii) using soft spins as well as hard 
spins. However, the corrections to scaling for small sizes 
always had the same sign. Consequently, we have not 
found a model where we could set the leading correction 
to zero by fine tuning a parameter in the Hamiltonian. 
We are not claiming that such a model does not exist; 
only that we were not able to find it. 

Second, include corrections to FSS in the analy- 
sis. Correction terms are characterized by an expo- 
nent, which is universal, and an amplitude which is not. 
Corrections to scaling in a scaling plot, where one at- 
tempts to collapse data from different sizes onto a single 
curve, occur for both the horizontal and vertical axes, see 
e.g., Ref. I30I Thus, a large number of additional param- 
eters have to be determined from the data when correc- 
tions are included. For ferromagnets, where extremely 
precise data can be obtained for a very large range of 
sizes, this is possiblei^ However, for spin glasses, the 
range of system sizes is more limited because of slow dy- 
namics, even though we have used state-of-the-art algo- 
rithms and considerable computer time, and the statis- 
tics are not as good because there are large variations be- 
tween different realizations of the disorder. Hence our at- 
tempted fits which included corrections to scaling did not 
determine the parameters well and frequently the nomi- 



TABLE I: Selection of different estimates (chronologically 
sorted) of the critical temperature Tc as well as the critical 
exponents computed by different groups for Gaussian (G) as 
well as bimodal (±J) random bonds. The estimates show 
strong variations and often do not agree. Note that Tc is not 
universal, so the issue at hand is whether or not the results for 
V and for 77 agree within the error bars. The results by Jorg 
(Ref. 113) are for a bond-diluted ± J spin glass with 45% bond 
occupation, which is why the estimate of Tc is different than 
for the standard bimodal spin glass. The results of Toldin et 
al. (Ref. [13) are for a random-anisotropy Heisenberg model 
(RA) in the strong-anisotropy limit, which is expected to be 
in the same universality class as the Ising spin glass in three 
dimensions. The last two rows describe results of the present 
study and will be described in detail in what follows. 



Authors 




Tc 


V 


n 


Ogielski & Morgensterr>~ 


±J 


1.20(5) 


1.2(1) 




Ogielskiii 


±J 


1.175(25) 


1.3(1) 


-0.22(5) 


McMiUaiii^ 


G 


1.0(2) 


1.8(5) 




Singh & Chakravarty— 


±J 


1.2(1) 


1.3(2) 




Bray & Mooreii 


G 


1.2(1) 


3.3(6) 




Bhatt & YoungiS 


±J 


1.2(2) 


1.3(3) 


-0.3(2) 


Bhatt & Youngia 


G 


0.95(5) 


1.6(4) 


-0.4(2) 


Kawashima & Young— 


±J 


1.11(4) 


1.7(3) 


-0.35(5) 


Bernardi et aim 


±J 


1.165(10) 




-0.245(20) 




G 


0.88(5) 




-0.50(4) 


Inigues et aim— 


G 


1.02(5) 


1.5(3) 




Berg & Janke— 


±J 


1.12(1) 




-0.37(4) 


Marinari et alw^ 


G 


0.95(4) 


2.00(15) 


-0.36(6) 


Palassini & Caracciolo«- 


±J 


1.156(15) 


1.8(2) 


-0.26(4) 


Mari & Campbell^ 


±J 


1.20(1) 




-0.21(2) 




G 


0.86(2) 




-0.51(2) 


Ballesteros et alw~ 


±J 


1.138(10) 


2.15(15) 


-0.337(15) 


Mari & Campbell^ 


±J 


1.190(15) 




-0.20(2) 




G 


0.920(15) 




-0.42(2) 


Mari & CampbellSa 


±J 


1.195(15) 


1.35(10) 


-0.225(25) 


Nakamura et alm~ 


±J 


1.17(4) 


1.5(3) 


-0.4(1) 


Pleimling & Campbell" 


±J 


1.19(1) 




-0.22(2) 




G 


0.92(1) 




-0.42(2) 


Jorgia 


±J 


0.663(6) 


2.22(15) 


-0.349(18) 


Campbell et alw~ 


±J 




2.72(8) 


-0.40(4) 


Toldin et alP 


RA 


0.93(4) 


2.4(6) 


-0.24(4) 


This study 


G 


0.951(9) 


2.44(9) 


-0.37(5) 




±J 


1.120(4) 


2.39(5) 


-0.395(17) 



nal "best" fit had extremely large corrections to scaling 
and unphysical values for the parameters. 

Since attempts to eliminate (leading) corrections to 
scaling and to explicitly incorporate them failed, we re- 
sorted to the strategy of just including data for the larger 
sizes where we expect corrections to be small, and ne- 
glecting corrections to scaling. Our main conclusion is 
that in equilibrium universality is satisfied, since results 
for a particular observable do not appear to depend on 
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the disorder distribution. However, there are still some 
open questions since different observables for a single dis- 
order distribution yield estimates of critical exponents 
which differ by more than the estimated (statistical) er- 
ror bars. This discrepancy can probably only be resolved 
by an analysis incorporating corrections to scaling. Al- 
though this does not appear to be possible at present 
(since the range of sizes is too small) it may be possible 
in the future if more significantly efficient algorithms can 
be developed. 

The paper is structured as follows: In Sec. we in- 
troduce the model as well as the measured observables 
and in Sec. IIIII we present details on the Monte Carlo 
method used. In Sec. |^ we discuss the method used to 
estimate unbiased error bars for the different estimates of 
the critical parameters. Results are presented in Sec. Ivl 
followed by concluding remarks in Sec. IVII 



II. MODEL AND OBSERVABLES 

A. Edwards-Anderson Model 

The Hamiltonian of the Edwards- Anderson Ising spin 
glassi*^ is given by 



n 



I J Si Sj , 



(1) 



where the sites i lie on a three-dimensional cubic lattice 
of size N = and the spins Si can take values ±1. The 
sum is over nearest-neighbor pairs and the interactions 
Jij are independent random variables. Periodic bound- 
ary conditions are applied. In this work we mainly study 
two paradigmatic cases of the EA model: 



been unsuccessful; therefore we have collected less good 
statistics for these models than for the Gaussian and ± J 
distributions. Hence we shall not present results for these 
other models in detail, except in Sec. IV Dl where we will 
do a global comparison of all the models studied to test 
for universality. These other models are 



• Gaussian/bimodal distribution with j) '^u ~ 0- 
This model is the same as the model presented in 
Eqs. Q and ©, but with the constraint that the 
sum of the Jij is exactly zero. 

• Correlated bonds: In this model the nearest- 
neighbor bonds have correlations. The probability 
distribution for the bonds is taken to be 



V (Jij) oc cxp 



77 Jij — ^^y^JiiJi2JaJn 



(4) 



where the last term involves the product of the four 
bonds around an elementary plaquette of the lat- 
tice, and is summed over all plaquettes. It is this 
term which generates correlations in the bonds. We 
take Jij = ±1, for which the first term in Eq. Q is 
actually a constant. We generate correlated bonds 
by first performing a Monte Carlo simulation on the 
bonds using the Hamiltonian in Eq. Q . The bonds 
are then frozen and the spin simulation is carried 
out. 

• Cosine disorder distribution: The bonds of the EA 
spin glass are chosen according to 



cos Jy 



(5) 



• Gaussian-distributed random bonds 
mean and standard deviation unity. 



with 



r{j,, 



1 



-J?,/2 



(2) 



• Bimodal (± J) distribution of bonds in which Jy 
take the values ±1 with equal probability. 

V{J,,) = ^[S{J,,-l) + SiJ,,+l)] . (3) 

In all cases that we study the mean of the distribution 
is zero. Since we set the standard deviation to be unity, 
the temperature is a dimensionless quantity. 



• Soft spins: In this model the spins Si can take any 
length from — oo to +oo and the Hamiltonian is 
given by^^ 



n 



ijSiSj 



Ysf + xYiSf-^", (6) 



where the bonds Jij are Gaussian distributed and 
the parameter A controls the average length of the 
spins. For A ^ oo we recover the Ising model with 
fixed-length spins. 



C. Measured Quantities 



B. Other Models 

We have also studied other models, although in less 
detail than the Gaussian and ± J models, in order to see 
if, by tuning parameters, we could eliminate the leading 
correction to scaling in any of them. These attempts have 



We measure different observables which, in the past, 
have proven to show a good signature of the phase tran- 
sition. First, we study the Binder cumulant'^^ given by 



4 



where (■ • • )t represents a thermal average, [••■]av is a 
disorder average, and q is the spin overlap given by 



Q ■ 



1 ^ 

4=1 



(8) 



In the previous equation "a" and "/?" represent two 
copies of the system with the same disorder. The Binder 
ratio is dimensionless and thus has the simple scaling 
form 



(9) 



where (3 — 1/T and /3c is the inverse of the critical tem- 
perature. In addition to /3c, the "metric factor"i^ A is 
also nonuniversal, but, since A is included explicitly, the 
resulting scaling function G(x) is universalMt For the 
models studied here with no lattice anisotropy^^^iSi^^ a 
universality class for finite-size scaling functions is spec- 
ified by (i) the hulk universality class, (ii) the boundary 
conditions, and (iii) the sample shape. In this work we 
always use the same boundary conditions (periodic) and 
the same sample shape (cubic), so we expect the same 
function G{x) for models which lie in the same bulk uni- 
versality class. According to Eq. data for g for differ- 
ent sizes should intersect at T^. Furthermore, the value 
of g at the intersection point, which is given by G'(O) is 
also universal since the whole function G{x) is universal. 
In addition, we study the finite-size correlation length 

^^24,25,36,37 ^ggj^g^ ^y 

-,1/2 



1 



2sin(fcmin/2) 



Xsg(O) 

XSG(kmin) 



1 



(10) 



where k,„in = (27r/L, 0, 0) is the smallest nonzero wave 
vector. Here, the wave vector dependent spin-glass sus- 
ceptibility is given by 



i,3 



(11) 



= X 



(al^'''[P-P, 



(12) 



Like the Binder ratio, the finite-size correlation length 
divided by the system size is a dimensionless quantity 
and so scales as 

L 

in which the metric factor A is the same^ as in Eq. 
The reason the metric factors are the same is as fol- 
lows: By hypothesis, the argument of all ESS functions 
is really i/^oo: with no metric factor, where is the 
bulk correlation length. In this form one has separate 
ESS functions for each side of the transition because 
l/^oo = B±\[3 — PcY" vanishes in a singular manner at 
criticality. Eor example, 

Y = ^± (i/^oo) (13) 
- X±{B±L[P-I3,r) (14) 

(15) 



in which X± and X± are universal. Comparing Eq. H15() 
with Eq. (O, we see that there are universal, i.e., dis- 
tribution independent^ factors c± such that 



(16) 



and therefore X and X are essentially the same functions, 
in the sense that 



X{u) = 



_ ] X+{u/c+) {u > 0), 



X_(u/c-) (ii<0). 



(17) 



If we repeat the argument for the Binder ratio we 
obtain 



(18) 



and choosing the same c± as in Eq. H16() we reproduce 
Eq. lO with the same value for A as in Eq. (|12|l . 

The advantage of £,l/L over the Binder ratio g is 
that the Binder ratio, being restricted to the interval 
g G [0, 1], does not have much room to "splay out" be- 
low Tc- Presumably because of this, the data for g in 
three dimensions depend only very weakly on the system 
size in this region i22i2^ However, the finite-size correla- 
tion length ^l/L does not have this constraint, and so 
the data for it splays out better at low temperatures, al- 
lowing for a more precise determination of the critical 
parameters. 

Both g and ^l/L allow one to determine Tc and the 
critical exponent v. However, to fully characterize the 
critical behavior of a system, a second critical exponent, 
1], is required.™ Thus we also study the scaling of the 
spin- glass susceptibility xsG = XSG(k = 0) given by 
Eq. I|ll|) with k = 0. Near criticality we expect 



XSG = DL^-^>ci^AL'/-[P- (3,]) 



(19) 



therefore allowing us to determine the critical exponent 
ry. By separating out the nonuniversal amplitude D, the 
scaling function C is universal. 



III. NUMERICAL DETAILS 

The simulations are done using the parallel tempering 
Monte Carlo methodj^S*^ Eor the Gaussian distribution 
we have tested equilibration with the method introduced 
m Ref. lill where the energy computed directly is com- 
pared to the energy computed from the link overlap. The 
data for both quantities approach a limiting value from 
opposite directions. Once they agree, and other observ- 
ables are independent of Monte Carlo steps, the system 
is in equilibrium. Eor the bimodal disorder distribution 
we use a multispin coded version of the algorithm which 
allows us to update 32 copies of the system at the same 
time. The aforementioned equilibration test cannot be 
applied to the bimodal spin glass. In this case we study 
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TABLE II: Parameters of the simulations for Gaussian dis- 
tributed disorder. N^^ is the number of samples, N^^ is the 
total number of Monte Carlo sweeps for each of the 2Nt repli- 
cas for a single sample, Tmin is the lowest temperature sim- 
ulated, and A'^r is the number of temperatures used in the 
parallel tempering method for each system size L. 



L 




Asw 


^min 


At 


3 


20000 


32768 


0.80 


8 


4 


20000 


20000 


0.80 


8 


6 


20000 


40000 


0.80 


8 


8 


20000 


50000 


0.80 


10 


12 


10000 


655360 


0.80 


16 


16 


5000 


1048576 


0.80 


33 


TABLE III: Parameters of the simulations, defined in Ta- 
ble ^ for bimodal distributed disorder. The system sizes 
marked with an asterisk have been simulated with the more 
efficient multispin method. 


L 


Afea 


Asw 


Tin ill 


At 


3 


40000 


8000 


0.82 


16 


4 


40000 


8000 


0.82 


16 


6 


40000 


20000 


0.82 


16 


8 


30000 


80000 


0.82 


16 


12 


15807 


300000 


0.82 


18 


16* 


11360 


128000 


0.95 


16 


20* 


9408 


1280000 


1.05 


25 


24* 


8416 


1280000 


1.05 


25 



how the results vary when the simulation time is suc- 
cessively increased by factors of 2 (logarithmic binning). 
We require that the last three results for all observables 
agree within error bars. 

Parameters of the simulation are presented in Tables UTI 
and mil for the Gaussian and bimodal (± J) distributions, 
respectively. 

IV. STATISTICAL ANALYSIS OF THE DATA 

For the Binder ratio and finite-size correlation length 
we need to find the best choice of parameters, v and /?c in 
order to collapse the data onto the scaling predictions of 
Eqs. © and ((T^ . respectively. To do so we assume that 
the scaling function can be represented by a third-order 
polynomial y{x) = cq + cix + C2X^ + c^x^ and do a global 
fit to the six parameters Ci, (i = 0, • • • , 3), /3c, and v. We 
also analyze results for xsG , for which there is a seventh 
parameter, the critical exponent r]. We have performed 
these nonlinear fits in two ways: (i) a code based on the 
Levenberg-Marquardt algorithm;^ and (ii) the statistics 
package R.-^ The same results have been obtained from 
both approaches. 

It is also necessary to obtain error bars on the fit pa- 
rameters. One has to be careful because, for a given 
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FIG. 2: (Color online) Finite-size correlation length ^l/L as 
a function of inverse temperature f3 for the three-dimensional 
Edwards- Anderson Ising spin glass with Gaussian disorder 
for several system sizes L. The data cross at P^T^ ~ 0.951. 
The dashed lines represent the optimal values obtained from 
a finite-size scaling for /3c and ^l{Pc)/L. 

size, all temperatures are simulated with the same dis- 
order realization in the parallel tempering Monte Carlo 
method. Hence the fitted data are correlated. We 
therefore have applied the following procedure: For each 
system size L with Nsa, disorder realizations, a randomly 
selected bootstrap^'* sample of A^sa disorder realizations 
is generated. With this random sample, an estimate of 
the different observables (with bootstrap error bars) is 
computed for each temperature. We repeat this proce- 
dure A^boot = 1000 times for each lattice size and then 
assemble A'^boot complete data sets (each having results 
for every size) by combining the «-th bootstrap sample for 
each size for i = 1, • • • , A^boot- The finite-size scaling fit 
described above is then carried out on each of these A'boot 
sets, thus obtaining A'boot estimates of the fit parameters. 
Since the bootstrap sampling is done with respect to the 
disorder realizations which are statistically independent 
we can use a conventional bootstrap analysis to estimate 
statistical error bars on the fit parameters. These are 
equal to the standard deviation among the A'boot boot- 
strap estimates. 



V. RESULTS 

A. Gaussian-distributed random bonds 

Figure |2] shows data for the finite-size correlation 
length as a function of the inverse temperature for dif- 
ferent system sizes. The data for L > 6 intersect at (or 
very close to) a common point whereas the data for the 
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FIG. 3: (Color online) Finite-size scaling analysis of ^l/L 
for the Ising spin glass with Gaussian bonds according to 
Eq. p2|l . The sizes are 6 < L < 16. The solid line is the 
third-order polynomial used in the fit, see Sec. liVI 



smallest sizes, L = 3 and 4 lie consistently too low in this 
region. The fact that sizes L = 3 and 4 do not intersect 
at a common point clearly indicates that corrections to 
scaling are significant for these sizes. We were hoping 
to find other models where the trend would be the other 
way around (i.e., where the small-L data are too high) so 
that by fine tuning of parameters we could eliminate this 
correction to scaling. However, all the models studied 
(see Sec. mi had corrections of the same sign as shown in 
Fig.H 

In scaling the data according to Eq. (|12|l in the way 
discussed in Sec. II VI we omit sizes L — Z and 4 because 
these data are clearly affected by corrections to scaling, 
and Fig. 121 shows the resulting plot for sizes L > 6. The 
overall fit is very satisfactory and gives the critical pa- 
rameters shown in Table HVl 

In Fig. 01 we show data for the Binder ratio, Eq. {Tj), as 
a function of the inverse temperature /3. The data cross 
at ~ 0.931. Note that for (3 > Pc the data do not 
splay very well thus making it difficult to determine the 
critical temperature accurately. 

Using the analysis presented in Sec. IIVI we have ob- 
tained the best fit, shown in Table Hvl and present the 
scaling plot of the Binder ratio in Fig[Sl For the analysis 
we have only considered L > 6. 

Overall, we expect that the analysis of ^l/L gives more 
accurate results than that for g, because the data for g 
do not splay out much below Tc, and so our best results 
for Tc and ly are those for ^l/L, i.e. 



Tc = 0.951(9), 



2.44(9) (Gaussian) . 



(20) 



TABLE IV: Summary of critical parameters for a Gaussian 
disorder distribution estimated by scaling the data. Scaling 
has been done in (/3 — /3c) except for the spin-glass suscepti- 
bility for which the data has been scaled with (T — Tc). In the 
table below, Tc is the critical temperature, and r] and v are 
critical exponents. The quantity co is the zeroth-order coeffi- 
cient of the fitting polynomial and corresponds to the value of 
a given observable at criticality. represents the chi-squared 
value for the finite-size scaling fitting function (Ref. l42t) . For 
comparison the number of data points used in the fit is 25. 
For the fit using the scaling form of Ref. the value of Tc is 
fixed to be that obtained from ^l/L. 





Estimate 


Error 


Co 


0.6346 


0.0090 


Tc 


0.9508 


0.0089 


V 


2.4370 


0.0924 




11.7859 


10.2696 


9 


Estimate 


Error 


Co 


0.7600 


0.0068 


Tc 


0.9310 


0.0137 


u 


2.6761 


0.1662 




17.7245 


14.0111 


XSG 


Estimate 


Error 


Tc 


0.9489 


0.0264 


V 


1.4859 


0.0602 


V 


-0.3733 


0.0483 


x' 


12.8776 


8.0025 


XSG (scaling as in Ref. 


Estimate 


Error 


Tc 


0.9508 


0.0089 


V 


2.7767 


0.0249 


V 


-0.3716 


0.0055 


x' 


18.3403 


12.6776 



which are hard to estimate for the range of sizes that can 
be studied. We discuss systematic errors below in more 
detail, especially in Sees. IV CI and IVII It is gratifying 
that the results obtained from the analysis of g, namely, 
Tc = 0.931(17), ly = 2.67(17) are (just) consistent with 
these. In addition to exponents, the values of ^l/L and g 
at Tc are also expected to be universal. These are given 
by the appropriate values of cq in Table Hvl 



We emphasize that the error bars quoted in this paper 
are only statistical. There are also systematic errors. 



^^^^ = 0.635(9), 5(Tc) = 0.760(70) (Gaussian) . 

(21) 

Unfortunately, the situation for the analysis of the 
spin-glass susceptibility data is less gratifying. The scal- 
ing plot for the spin-glass susceptibility is shown in Fig.El 
For consistency with the other plots the horizontal axis 
in this figure is p. However, the method of fitting (third- 
order polynomial) works best, in this case, for a fit using 
T. Hence Fig. |H1 indicates the fit parameters from the T 
fit. 
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)S 

FIG. 4: (Color online) Binder ratio g, defined in Eq. JTJ, as 
a function of inverse temperature P for the three-dimensional 
Ed wards- Anderson Ising spin glass with Gaussian bonds. The 
data cross at (3^^ ~ 0.931. The dashed lines represent the 
optimal values obtained from a finite-size scaling for /3c and 



-1.5 -1 -0.5 0.5 1 1.5 

FIG. 6; (Color online) Finite-size scaling of the spin-glass sus- 
ceptibility XSG according to Eq. I19I I for Gaussian disorder. 
In fact, the method of scaling worked better when using T 
rather than /3 and the fit parameters shown are for the T fit. 
However, we show the resulting scaling plot using (3 for consis- 
tency with the other plots. The scaling analysis is performed 
for L > 6. 
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FIG. 5: (Color online) Finite-size scaling analysis of the data 
for the Binder ratio g according to Eq. for the three- 
dimensional Ising spin glass with Gaussian disorder. The 
scaling analysis is performed for L > 6 and the solid line 
represents the best fit to the data from the finite-size scaling 
analysis. 

The results of the fit are 

Tc = 0.949(26), 1^=1.49(6), ?? = -0.37(5). (22) 

The value for Tc agrees within the error bars with those 
from ^l/L and g, but the value for v is in strong dis- 



agreement. Disagreements between exponents obtained 
in different ways are presumably due to corrections to 
FSS, but the size of the difference here is surprisingly 
large. In Sec. IV CI we shall revisit the problem of the 
surprisingly low value for i/ obtained from xsG- 



B. Bimodal-distributed random bonds 

We show data for ^l/L in Fig.[7| It is fairly similar to 
the data for the Gaussian distribution in that the larger 
sizes show a common intersection, but the data for the 
smaller sizes, L = 3 and 4, are lower, showing that these 
sizes are affected by corrections to FSS. Hence we only 
use sizes 6 < L < 24 in the scaling plot which is shown 
in Fig. |S1 Parameters obtained from the fits are shown 
in Table El 

Figure El shows data for the Binder ratio, Eq. O, for 
different system sizes i as a function of the inverse tem- 
perature p. The corresponding scaling plot is shown in 

Fig.nni 

The fit parameters obtained from ^l/ L are 

Tc = 1.120(4), 1/ = 2.39(5) (± J) . (23) 

Those obtained from g, namely, Tc = 1.088(6), v — 
2.79(11) disagree somewhat, but, as also discussed above, 
we feel that those for (,l/L are more reliable because the 
data for £,l/L splay out more below Tc. 

The values of Cl/L and g at Tc are given by the ap- 
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TABLE V: Summary of critical parameters for a bimodal 
disorder distribution estimated by scaling the data. The scal- 
ing is done in (/3 — /3c) except for the data for the spin-glass 
susceptibility where the scaling is done in (T — Tc). For fur- 
ther details see the caption of Table Hvl The number of data 
points used in the fits is 48. In the fit for XSG using the scal- 
ing form of Ref.^^the value of Tc is fixed to be that obtained 
from ^l/L. 



Estimate 



Error 



Co 
1/ 



0.6265 
1.1199 
2.3900 
52.8369 



0.0036 
0.0037 
0.0514 
28.3532 



Estimate 



Error 



Co 
1/ 



0.7626 
1.0881 
2.7937 
60.5020 



0.0029 
0.0062 
0.1103 
30.2953 



XSG 



Estimate 



Error 



V 

V 



1.1040 
1.5721 
-0.3954 
88.2526 



0.0097 
0.0251 
0.0168 
31.8466 



XSG (scaling as in Ref. Estimate Error 

V 



1.1199 
2.7376 
-0.3663 
83.6070 



0.0037 
0.0166 
0.0166 
36.7894 
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FIG. 7: (Color online) Finite-size correlation length £,l/L as 
a function of inverse temperature f3 for the three-dimensional 
Ed wards- Anderson Ising spin glass with ± J bonds for several 
system sizes L. The data cross at p-^ ~ 1.12. The dashed 
lines represent the optimal values obtained from a finite-size 
scaling for /3c and ^l{Pc)/L. 
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FIG. 8; (Color online) Finite-size scaling analysis oi ^l/L for 
the Ising spin glass with ±J bonds according to Eq. I12II . The 
scaling analysis is performed for L > 6. 
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FIG. 9: (Color online) Binder ratio g as a function of inverse 
temperature /3 for the three-dimensional Edwards- Anderson 
Ising spin glass with ±J bonds for several system sizes L. 
The data cross at f]^^ — 1.088. The dashed lines represent 
the optimal values obtained from a finite-size scaling for /3c 
and g{Pc)/L. 



propriate values of cq in Table IVl 

^I^^ = 0.627(4), g{Tc) = 0.763(3) (± J) . (24) 

A scaling plot for xsG is shown in Fig. ^] The best 
fit (using T rather than (3) gives 



1.104(9), j/ = 1.57(3), = -0.395(17). (25) 
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FIG. 10: (Color online) Finite-size scaling analysis of the data 
for the Binder ratio g according to Eq. Q for the three- 
dimensional Ising spin glass with ±J bonds. The sizes used 
are 6 < L < 24. 




FIG. 11: (Color online) Finite-size scaling of the spin- 
glass susceptibility XSG according to Eq. 1191 1 for the three- 
dimensional Ising spin glass with ±J bonds. As for the cor- 
responding plot for the Gaussian distribution, Fig. |S1 the fit 
is actually done in T, but the plot is given as a function of /? 
for consistency with the other plots. The scaling analysis is 
performed for L > 6. 



The value of differs from that obtained from ^l/L, in 
Eq. (|23|l by rather more than the error bars. If we fix Tc 
to be that obtained from ^l/L we find 

V = 1.527(8), 7? = -0.368(24). (26) 

The value of rj obtained in this way is consistent with 



that in the unconstrained fit in Eq. (|25|l . The value of 
is slightly different from that in Eq. H25|) . but more 
importantly, both these values of v obtained from XSG 
are considerably smaller than those obtained from /L 
and g. This is the same situation than found for the 
Gaussian distribution. Interestingly, the values of i> from 
XSG for the two distributions agree quite well with each 
other. 

We argue that the most reliable quantity to analyze 
is £,l/L because this has clean intersections with signifi- 
cant splaying out below Tc. The results for the exponent 
I/, shown in Eqs. H2U|) and (|23|l . for the Gaussian and 
± J distributions agree well within the error bars, which 
supports universality. Further support for universality 
comes from the agreement in the values of / L and g at 
the critical point, shown in Eqs. (|21|) and H24|) . We also 
note the agreement in the values of rj from Eqs. (|22|l and 
(ESJ) or P^. 



C. Alternative analysis of Xsg 

The main unresolved issue is the large difference in 
the values of v obtained from the spin-glass suscepti- 
bility compared with those obtained from S^l/L and g. 
Recently, Campbell et alm^ have claimed that the differ- 
ence is much diminished if one uses an alternative scaling 
form. They propose that the scaling region will be larger, 
so one can incorporate data for a larger range of temper- 
ature, if the behavior as T ^ oo is consistent with the 
scaling function. To be precise, they propose that 
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G 



{LTY''' (1 - {TjTf)\ , (27) 
'{LTf- (1 - (Te/r)2)] , (28) 
XSG = {LTf-^c\{LTf/'' {l-{TjTf)\, (29) 



L 



where we have not included explicitly the metric factors. 
Asymptotically, for L — > oo and (T — Tc) ^ 0, these ex- 
pressions are equivalent to the standard forms that we 
have used, Eqs. ©, ifT^ . and 1(1^)1. Thus the difference 
between the expressions proposed by Campbell et al. and 
the standard expressions is only in the corrections to scal- 
ing. 

In both the original scaling forms and the modified 
form of Campbell et al., Tc is located by intersections 
of data for ^l/L and g of different sizes. Thus we do 
not expect the estimates of Tc to be very different in 
the two approaches. Furthermore, if we restrict data to 
the region close to Tc, the data collapse involves mainly 
the derivative of the data with respect to temperature at 
Tc. For ^l/T and g, both the original and modified form 
predict that the temperature derivative is proportional 
to L^/'^ . Hence, for ^l/T and g, we also do not expect 
very different values for v from the two scaling forms. 
These expectations are confirmed by our analysis. Using 
Eqs. H27|) and H28|l . we find values for Tc and v which 
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agree, within the error bars, with those described above 
in Sees. IV Al and IV ij| which used the standard scahng 
forms. It is possible that scahng may work for data over 
a larger range of (T — Tc), at least above^^ Tc, using 
Eqs. (|?7jl and but we have not investigated this in 
detail. 

However, the situation for xsG is quite different, be- 
cause of the factor of T'^~^ in front of the scaling function 
in Eq. since 



1 dxsG 



XSG dT 



(30) 



where a and b depend on Tc and the value of the scal- 
ing function and its derivative at zero argument, but not 
on L. The factor of b arises from the T dependence of 
the prefactor outside the scaling function in Eq. (|29|l , and 
does not occur in the analogous expressions for L and 
g in Eqs. |(2H1) and (|2I|). For L ^ oo, the L^/'' term in 
Eq. H3U|I dominates and we recover the same behavior 
as in the standard scaling form. However, for the small 
sizes that can be studied numerically, the factor of b gives 
a significant correction to scaling especially since l/v is 
small. This is why Campbell et alm^ found a large differ- 
ence in the value of v obtained from from xsG using their 
scaling compared with conventional finite-size scaling. 

We also find a large difference. If we fix to be the 
value obtained from £,l/L we obtain, from Eq. H29|) 



J/ = 2.777(25), T] 
V = 2.738(17), 77 



-0.372(6) 
-0.366(3) 



(Gaussian), (31) 
(±J), (32) 



see also Tables HVI andlVI These values for v are consider- 
ably larger than those found from the standard analysis 
discussed in Sees. IV Al and IV Bl They are even some- 
what larger than those found from S^l/L, although they 
agree better with the £,l/L values than those found from 
XSG using the standard analysis. The fact that we, like 
Ref. |2^ obtain very different values for v from xsG de- 
pending on the form of the scaling function used, tells 
us that corrections to scaling can be very important in 
spin glasses for the range of sizes that can be simulated. 
We note, however, that the two estimates in Eqs. H31|) 
and (|32() agree well with each other, so we still find no 
evidence for lack of universality. 



D. Global comparisons of all the models 

We have computed two dimensionless quantities, L 
and which intersect at a finite value at the critical 
temperature. In the previous parts of this section, we 
have plotted both of them against (3. It turns out also 
to be useful to plot one of them against the otherm^^^ 
According to Eq. AL^I^IP - /3c] = X-^^/i) and 
so, from Eq. 0, we can write 
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FIG. 12: (Color online) Binder ratio g as a function of the 
finite-size correlation length L for several system sizes L > 
6 and for all models described in Sec. |n] All data collapse 
onto a universal curve, thus providing clear evidence that spin 
glasses obey universality. Different system sizes are labeled 
with different colors, and different models use different types 
of symbols as indicated. 



where G is a universal function. Note that there are 
no nonuniversal metric factors in this expression. Hence 
data for all the models described in Sec.llllshould collapse 
when g is plotted against £,l/L. This works very well as 
shown in Fig. 1121 which includes sizes L>Q. Figure IT^ 
provides additional very strong evidence for universality 
in spin glasses. 



VI. SUMMARY AND CONCLUSIONS 

We have studied numerically the phase transition in 
a variety of Ising spin-glass models in three dimensions, 
to test for universality. Our most detailed simulations 
are on nearest-neighbor models with Gaussian and ±J 
interactions, and our results for them are summarized in 
Tables IIVI and A comparison shows that correspond- 
ing estimates for the exponents v and 77 agree well, as 
do the values of ^l/L and g at criticality (labeled cq). 
This supports universality, as does the plot of g against 
i,L/L, Fig. El where data for all the models studied (not 
just the Gaussian and ± J models) collapse onto a single 
universal curve. 

The main unresolved issue is the large difference be- 
tween the values for v obtained from £,l/L or g on the 
one hand and xsG on the other. This is presumably due 
to systematic errors coming from corrections to scaling, 
but unfortunately we have not been able to incorporate 
corrections in our analysis since we do not have data with 
sufficient precision over a sufficiently large range of sizes. 
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The errors quoted in this paper are statistical errors only; 
systematic errors are not included. Evidence for strong 
corrections was found explicitly in Sec. IV CI where we 
used a scaling form for XSG proposed in Ref. 28 which 
differs from the standard form only in corrections to scal- 
ing. From this scaling form, we obtain an estimate for i' 
from our data for xsG which is very different from that 
obtained from xsG using the standard analysis. This 
large difference in the values of v from the two methods 
of analysis does not occur, however, for our data for g 
and ^l/L. 

Overall, we have found no evidence for lack of uni- 
versality, but have found evidence for strong correc- 
tions to s caling. We suspect that the dynamical data 
of Refs. which was interpreted to show lack of 

universality, more likely shows evidence for corrections 



to scaling. 
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